Sampling the ground-state magnetization of d-dimensional p-body Ising models 
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We demonstrate that a recently introduced heuristic optimization algorithm [Phys. Rev. E 83, 
046709 (2011)] that combines a local search with triadic crossover genetic updates is capable of 
sampling nearly uniformly among ground-state configurations in spin-glass-like Hamiltonians with 
p-spin interactions in d space dimensions that have highly degenerate ground states. Using this 
algorithm we probe the zero-temperature ferromagnet to spin-glass transition point q c of two example 
models, the disordered version of the two-dimensional three-spin Baxter-Wu model [q c = 0.1072(1)] 
and the three-dimensional Edwards- Anderson model [q c = 0.2253(7)], by computing the Binder 
ratio of the ground-state magnetization. 

PACS numbers: 75.50.Lk, 75.40.Mg, 05.50.+q, 64.60.-i 



I. INTRODUCTION 

Disordered systems with degenerate ground states of- 
ten have an exponentially-largei number of states that 
minimize the Hamiltonian (cost function) of the problem. 
A hallmark model system to study degenerate ground 
states is the Edwards Anderson Ising spin-glass model 
with bimodal (± J) interactions^ for which the entropy 
is extensive even at zero temperature due to the expo- 
nential number of ground-state configurations. 

Studying these systems at finite temperatures using 
Monte Carlo simulations typically poses no major chal- 
lenges: An estimate of an observable (e.g., the magnetiza- 
tion) has to be measured and averaged over Monte Carlo 
time. This means that an average over different degen- 
erate states is automatically taken into account. How- 
ever, when studying the ground-state properties of such 
systems, it is imperative to ensure that an average over 
either all ground-state configurations or at least an unbi- 
ased subset is taken. While some algorithms are known 
to do a relatively fair sampling of the ground-state con- 
figurations (e.g., simulated annealing),—^ others, such as 
quantum annealing^— have been shown to favor certain 
configurations and exponentially suppress others. 

Although finding ground-state energies is not typi- 
cally harder with discrete energy distributions than in 
models with continuous energies, determining thermody- 
namic quantities requires thorough checks of the algo- 
rithms used. Numerically exact results are achievable 
only in limited scenarios^— though comparisons with 
theoretical predictions also may provide confidence that 
the technique is behaving correctl y 12 ' 13 

Here, we investigate the effectiveness of a heuristic 
optimization algorithm previously developed^ for find- 
ing all ground states of a system, or, if that number is 
too large, finding a representative and unbiased sample, 
applied to a p-spin generalization of the ti-dimensional 
Edwards- Anderson Ising spin-glass modeL 15 ' 16 The gen- 
eral applicability of this algorithm implies that it is use- 
ful for studying a wide range of hard problems where 
it is desirable to average over many degenerate opti- 



mal cases. We illustrate the usefulness of the method 
for studying spin-glass models by computing the zero- 
temperature fcrromagnet-to-spin-glass phase transition 
as the disorder strength is varied in the disordered 
two-dimensional three-spin Baxter-Wu model and in 
the three-dimensional two-spin Edwards- Anderson model 
with bimodal interactions. 

Using a moderate numerical effort, we show for the case 
of the three-dimensional two-spin Edwards- Anderson 
Ising spin-glass model that the critical threshold q c be- 
tween a ferromagnetic and a spin-glass phase is q c = 
0.2253(7), in agreement with previous studies^ albeit 
with considerably smaller error bars. 

In Sec. |TT] we present an overview of the algorithm, fol- 
lowed by results on the two-dimensional three-spin model 
in Sec. lIIII and results on the three-dimensional Edwards- 
Anderson spin glass in Sec. |TV] 



II. NUMERICAL METHOD 

The Edwards Anderson model in three space dimen- 
sions has been shown to be NP-hardr^ as has the spin 
glass with three-body interactions in two dimensions^ It 
is therefore expected that there are no exact algorithms 
capable of simulating large instances of either model at 
zero temperature. The genetic algorithm presented in 
Ref. [Til heuristically generates many optimal ground- 
state configurations of disordered p-spin models. The 
method starts from a population of N p random spin con- 
figurations and combines a highly-effective local search 
heuristic with triadic crossover genetic updates^ Al- 
though very large systems are unattainable with such 
an approach, it was shown to produce ground-state en- 
ergies in the two-dimensional three-spin model with high 
confidence for N < 42 2 = 1764 spins-ii 

Genetic algorithm: The algorithm proceeds according 
to the following steps 

1. Initialize N p configurations: for each configuration 
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(a) Initialize the spins to a high-temperature con- 
figuration where all spins are assigned ran- 
domly. 

(b) Perform a local search update to generate a 
low-energy configuration. 

(c) If this configuration is already present in the 
population, go back to 1(a). 

(d) Otherwise, accept this configuration and pro- 
ceed to the next one. 

2. Perform a triadic crossover update to generate two 
new configurations, C± and C 2 . 

(a) Select three parents Pi, P 2 and P3 randomly. 

(b) For each spin i 

i. If Pi(i) = P 3 {i), Ci(i) <- Pi(i) and 
C 2 (*)^P 2 (z) 

ii. Otherwise Ci(i) P 2 (i) and C 2 (i) 
Mi) 

(c) Select two configurations X and 1" from the 
population at random. 

(d) If E(d) < E(X), replace X with d; if 
P(C 2 ) < E{Y) replace F with C 2 . 

3. Repeat step 2 iV s teps times or until all states of the 
system have the same energy. 

Unlike in many genetic algorithms, mutations are not 
necessary to achieve good results with this method. 

Local search component: The local search algorithm 
used in step 1(b) above (described in detail in Ref. [Til ) 
generates clusters of nearby spins for which the energy 
is lowered if all spins in the cluster are flipped. This is 
achieved by performing a depth-first search starting from 
each spin in the system. The search algorithm takes two 
parameters, E max > and d max . At each step of the 
search, the spins adjacent to the current site that are not 
already in the cluster are considered to be added. As each 
spin is added, the energy of flipping the current spin and 
all its ancestors in the search tree (all current members of 
the cluster) is maintained. If the energy is negative, the 
algorithm has found an energy-lowering move, so all spins 
in the current cluster are flipped and the search termi- 
nates successfully. If the energy exceeds the threshold en- 
ergy -E max or the search depth becomes more than d max , 
this search branch is discarded. For E max and d max both 
large, this search algorithm gives strong results, but runs 
slowly. For use as part of the genetic algorithm above, 
speed is more important than finding the lowest energy, 
therefore we search for small clusters with few barriers, 
setting E max = 4 (the smallest possible change in energy 
in our model) and d max = L, where L is the linear size 
of the system. 

With N p and N stC p S sufficiently large, this genetic algo- 
rithm produces ground states with high confidence.— In 
general, to study a ferromagnet-paramagnet transition, 



it is useful to study the finite-size scaling of dimension- 
less quantities based on the magnetization of the system. 
By computing quantities such as the Binder ratio,— one 
can straightforwardly measure the location of the critical 
point, as well as the critical exponent v. When the disor- 
der distribution is continuous, each disorder sample has a 
unique ground state (up to global degeneracies), so find- 
ing this ground-state configuration immediately yields 
both the energy and magnetization; computing its energy 
and magnetization are therefore typically of comparable 
difficulty. But in the discrete-disorder case, the number 
of ground-state configurations N g is typically very large. 
Computing the magnetization requires either finding all 
ground states, or generating a uniform sampling to find 
a subset of these. This heuristic algorithm does not stop 
after it finds one ground-state configuration; it continues 
finding others until either all ground states are found or 
the algorithm can find no more ground states due to a 
built-in constraint. 

For a genetic algorithm, the quality of the results de- 
pends both on the run time of the algorithm and the 
number of configurations being concurrently updated, 
the population size N p . For sufficient run time, our 
numerical results show that this algorithm finds all the 
ground-state configurations of a degenerate system if the 
number of ground states N g < N p . When N g > N p , 
which will inevitably be the case for some large systems, 
the probabilities for finding each ground state vary some- 
what, although our tests suggest that the probability of 
finding a particular ground state is within some factor 
times the expected probability with a uniform distribu- 
tion, with that factor typically being of order unity. For 
one instance of the 3-spin model studied in Sec. IIIII with 
600 ground states, we reran our algorithm 10 5 times to 
find the relative probability for the genetic algorithm to 
find each ground state. Note that disorder instances with 
this number of ground states are rare for L = 6; the me- 
dian number of ground states we find over all disorder in- 
stances is 4. For this sample, when N p = 2N g all ground 
states were found in every run, and for N p = N g almost 
all ground states were found in nearly every runs. When 
N p = 0A8N g , the ground state with the largest basin of 
attraction is roughly eight times more likely to be seen 
than the ground state with the smallest. To probe the 
system size dependence, we also looked at a larger sys- 
tem of size L = 18, where we chose a sample that also 
has 600 ground states. In this case, the algorithm selects 
among those ground states much more uniformly, with 
only a factor of 2 between most and least likely ground 
states. For L = 18, the median number of ground states 
is 768; with N p > 5000, a large majority of the disor- 
der instances considered have N g < N p and are expected 
to be treated fairly by the algorithm. These results are 
shown in Fig. [T] 

Even when N g > N p , this performance strongly out- 
performs quantum annealing^— i2I where in the case of 
the fully-frustrated Villain models certain ground-state 
configurations are exponentially suppressed^ We have 
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FIG. 1: (Color online) Normalized relative probability p of 
each ground state occurring for varying population sizes. A 
single instance of disorder with N = 6 x 6 spins is considered, 
in addition to one instance of disorder with TV = 18 x 18. The 
disorder cases were selected, for comparison, to have the same 
number of ground states, with N g = 600. The ground states 
(indexed by i) are sorted according to their probability of oc- 
curring, giving a monotonically increasing plot. Even with 
N p < N g , the probability for finding each ground state is a 
factor of order unity times the uniform probability, implying 
that one can do a reasonable approximation of the T = 
thermodynamic state even when N p < N g . With other fac- 
tors being equal, the larger system size shows a more uniform 
sampling. 

not found any trends in the magnetizations produced by 
the algorithm; using only this magnetization data, the 
ground states appear to be chosen in random order, even 
when N g > N p , as is shown in Sec. lIIIl and Fig. [2] An av- 
erage over the configurations that have been found there- 
fore appears to be representative of the entire ensemble 
for this algorithm. 

III. THREE-SPIN ISING MODEL 
A. Model 

The three-spin Ising model with random plaquette in- 
teractions is given by the Hamiltonian 

'h = -jJ2vaS 1 a s 2 a s 3 a , (i) 

A 

where the sum is over all triangular plaquettes in a two- 
dimensional triangular lattice and each term involves the 
product of the Ising spins S A <G {±1} on the vertices of 
each plaquette. Without any loss of generality we set J = 
1. The plaquette couplings r] A are chosen independently 
and randomly according to a bimodal distribution, that 
is, 

P( VA ) = (1 - q)S(r ]A - 1) + q6( VA + 1) (2) 



with q — being the pure ferromagnetic case. As q in- 
creases, the model goes through a phase transition from 
ferromagnet to paramagnet. A previous study at finite 
temperature using Monte Carlo methods finds this tran- 
sition to occur at q c = 0.109(2)32 

This model is of interest for a number of reasons, p- 
spin models have been related to structural glasses^ - — 
Furthermore, when computing the error threshold of 
topological color code a 23 ' 28 to bit-flip errors the problem 
maps 2 ^ onto the 3-spin model in two space dimensions. 
Although the error threshold is given by the crossing 
of the ferromagnetic-paramagnetic phase boundary line 
with the Nishimori line^ the zero-temperature phase 
boundary delivers a lower bound for the threshold that 
is often easier to compute than to equilibrate finite- 
temperature Monte Carlo simulations. 

The phase diagram for the Ising model on a square 
lattice shares much in common with the three-spin Ising 
model on a triangular lattice. The three-spin Ising model 
with uniform ferromagnetic interactions is exactly solv- 
able and known as the Baxter- Wu model^ and the tran- 
sition temperature on a triangular lattice is identical to 
that of the two-dimensional Ising model on a square lat- 
tice. Although the models are in different universality 
classes; 23 ' 32 the two phase diagrams are quite similar even 
in the presence of disorder. Numerically, it appears that 
the multicritical points for these two models are very 
close to one another or identical £2- The results presented 
here suggest that this similarity does not extend to the 
low-temperature regime below the Nishimori line^ with 
possible reentrant behavior being much weaker, if it ex- 
ists at all, in this model. 

In contrast to the Ising model with p = 2, which has a 
two-fold (global spin flip) degeneracy, this model is four- 
fold degenerate. The magnetization 

M =^E^ (3) 

i 

in the ordered state is not symmetric about zero, al- 
though the average over the four degeneracy sectors re- 
lated by symmetry operations is zero: in the pure fer- 
romagnet, one of the four degenerate ground states has 
M = 1, while the other three states each have M = —1/3. 
In the paramagnetic state there is no state which is much 
stronger than the others, and the disorder-averaged mag- 
netization in each of these four sectors is zero. 

The zero-temperature thermodynamic state consists of 
all ground states with equal probability. To compute 
properties related to the magnetization, it is therefore 
not sufficient to find an arbitrary specific ground state of 
the model unless it can be demonstrated that it is chosen 
in an unbiased way. 

B. Sampling results 

When N p > N g , it is seen that the number of ground 
states found for the three-spin model is always a multiple 
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of 4, and the average magnetization is 0. It is necessary 
that these both be true due to the fourfold degeneracy 
of the system. This result suggests that the algorithm 
is finding all ground states of the system. It is possi- 
ble that there are other states which are not found, but 
we emphasize that the algorithm is given no information 
about the symmetries of the problem, so states which are 
related to one another by global spin-flip symmetries are 
not necessarily easy for the algorithm to link together. 
Thus is appears that we are seeing all ground states when 
N p > N g . 

When the number of ground states is large, so that 
N g > N p , it also appears that the subset of all ground 
states that the algorithm finds is representative of the en- 
tire ensemble. To test this, we study the magnetization 
of the ground states. The results for two typical cases 
are shown in Fig. [2] where the magnetization is shown 
for the first 5000 ground states found in the order i that 
they are found by the algorithm. There are no apparent 
temporal correlations in the output and the running av- 
erage quickly approaches zero, although there are some 
statistical fluctuations about this value because not all 
ground states have been found. This, in turn, justifies 
the use of this algorithm for sampling. 

The two samples shown in Fig. [2] are for two differ- 
ent instances of disorder with the same disorder strength 
q = 0.105 and the same system size (N = 18 2 ). This 
disorder strength is near the ferromagnet to spin-glass 
phase transition but slightly below it, so that the sys- 
tem is still ferromagnetic. In the sample shown in the 
top panel, the states all resemble the pure three-spin 
spin-glass problem, with 1 /4 of the configurations hav- 
ing M ~ 1, and the other 3/4 having M « —1/3. Unlike 
in the pure case, there are many distinct ground states, 
but each one follows the expectations for a ferromagnetic 
system. The sample depicted in the bottom panel shows 
the onset of disorder: while the majority of the states 
still have M ~ 1 or M » —1/3, a sizable fraction of the 
states here have M slightly below 0, and still others are 
between and 1/2. This implies that there are large col- 
lective excitations with zero-energy, which is a signature 
of glassy/frustrated systems. Thus, even for the same 
disorder strength, the system will behave ferromagncti- 
cally for some instances of the disorder and appear to be 
more glassy for others. 

It is instructive to compare against the two- 
dimensional Ising spin-glass model with p = 2. Even 
though efficient algorithms have long been known for 
computing ground states in this case^ these tech- 
niques have not been directly useful in finding all of the 
exponentially-many ground states, or even in finding un- 
biased randomly chosen ground states. Heuristic mea- 
sures have therefore been employed to find the ground- 
state magnetization^ although an efficient algorithm for 
fair sampling of the ground states for this model has been 
developed recently^ This model has been shown to ex- 
hibit reentrance from both zero- and finite-temperature 
simulations ; 33 ' 35 ' 36 and we are not aware of previous work 
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FIG. 2: (Color online) Left panels: Magnetization of 
ground-state configurations M for two independent samples 
(top/bottom panels) in the order they are found (indexed by 
i) for a three-spin model system with 18 2 spins and q = 0.105. 
The solid line is a running average that converges toward zero, 
suggesting that the ground-state magnetization is being sam- 
pled relatively fairly despite the slight bias shown in Fig. \T\ 
Right panels: Histograms P(M) of the magnetization. 



on the three-spin model below the Nishimori line for com- 
parison. 



C. Ferromagnet— to— spin-glass transition 

Although this algorithm does not permit the compu- 
tation of ground states with high confidence in very large 
systems, the intermediate system sizes accessible are ad- 
equate to precisely determine the location of the T = 
ferromagnet-to-spin-glass phase transition as the disor- 
der strength q varies. To do this, we compute the Binder 
rations for a restricted set of magnetizations. Due to the 
details of the four-fold degeneracy in this model including 
the lack of a global spin-flip symmetry, the ferromagnetic 
state consists of three states with weak negative magne- 
tization for every highly magnetized state with positive 



TABLE I: Summary 
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size N p and averaging 
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steps of q s tcp. 



of simulation parameters for the two- 
model of size L x x L y , with population 
; over N s different disorder instances, 
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magnetization. The paramagnetic state has no sector 
with stronger magnetization than the others. The fluctu- 
ations are strongest in the ferromagnetic sector, because 
there is less change from the behaviors in the other sec- 
tors between the two phases, so it is natural to focus on 
only the 1/4 of states that have a large magnetization. 
In Fig. [21 these are the states with M near 1. Let 



and 



M = M 



M k =M 



(4) 



(5) 



for numbered colors k = 1, 2, 3 of the tripartite lattice, 
with = 1 if site i has color k and —1 otherwise. We 
focus on the state of the system with the largest magne- 
tization, calling this the restricted magnetization 



M R = max(M ,M 1 ,M 2 ,M 3) 



(6) 



For any given spin configuration, it is straightforward to 
compute all four of these magnetizations. This is anal- 
ogous to taking the absolute value for the p = 2 Ising 
model, such that only configurations in the same state 
are compared against one another. We then compute 
the restricted Binder ratio <?r defined by 



9R 



Wk)U 

K M r)]Iv 



(7) 



where [■ ■ ■ ] av is an average over samples and (• ■ • ) rep- 
resents an average over ground-state magnetizations for 
the sample. This is the same as the standard definition of 
the Binder ratio except that the restricted magnetization 
Mr is used in place of M. For the p = 2 Ising model, 
taking the absolute value would not make any difference, 
but here the symmetric states have different magnetiza- 
tions. The restricted Binder ratio is dimensionlcss and is 
1 in the ferromagnetic phase and in the paramagnetic 



FIG. 3: (Color online) Scaling collapse of the restricted Binder 
ratio pa according to Eq. JS} for the two-dimensional 3-spin 
Ising model with system sizes L = 6 - 36. The dashed line is 
a curve fit to a third-order polynomial. Optimal data collapse 
is obtained for q c = 0.1072(1) and v = 1.5(1). 



phase (here, the T = spin glass). In a finite system of 
linear scale L with disorder strength in the vicinity of q Cl 
(7r is expected to scale as 



r M = G[L 1 ^(q-q c )]. 



(8) 



Equation (jSJ) allows for the extraction of the critical disor- 
der strength q c as well as the correlation-length exponent 
v. To estimate the optimal values of the critical param- 
eters we perform a finite-size scaling analysis where we 
tune q c and v until the chi-squared of a fit to a third- 
order polynomial in the vicinity of L 1 / U {q — q c ) < 1 is 
minimized. Error bars are determined by a bootstrap 
analysis. 37 ' 38 Our best estimates of the critical parame- 
ters are 



q c = 0.1072(1) 



1.5(1) 



(9) 



The scaling collapse using these critical parameters is 
shown in Fig. [3J To verify our results, we have also 
computed the standard Binder ratio on the unrestricted 
magnetizations from the same data, which produced con- 
sistent results for both q c and v. The only noticeable dif- 
ference in the results is the form of the scaling function, 
that is, gn(q = q c ) « 0.955, while g(q = q c ) w 0.23. A 
summary of the simulation parameters is shown in Ta- 
ble n 

Finally, to treat the possibility of finite-size corrections 
to scaling, we separate the components of the fit in Fig. [3] 
into L-2L pairs. We perform the curve-fitting procedure 
to extract q c using each pair of system sizes. The tri- 
angular lattice is tripartite only for L x mod 3 = and 
L y mod 2 = 0, so these comparisons are only possible for 
and N = LxL triangular lattice if L is a multiple of 6. To 
increase the number of cases we can handle, we consider 
L x to be any multiple of 3, and if L x is not divisible by 
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FIG. 4: (Color online) Measured values of q c as a function 
of system size. Each data point is the value of q c extracted 
from the scaling collapse of a system of size L with a sys- 
tem of size 2L, for L = 6, 9, 12, 15, and 18. The yellow 
(light) region is the estimate q c = 0.109(2) from Ref. [231 from 
finite-temperature Monte Carlo (MC) simulations of the 3- 
spin model. The solid green (darker) line is the T = esti- 
mate q c = 0.1072(1) in this study, and the blue (dark) region 
is, for comparison, the T = estimate q c — 0.103(1) for the 
two-dimensional 2-spin Ising model from Ref. [3^. The fact 
that the data (solid green line) and the result for the two- 
dimensional Ising model (blue/dark region) do not overlap 
clearly shows that while both models apparently share the 
same multicritical point, the behavior for temperatures below 
the multicritical point is quite different. 



2, we simulate for both L x x (L x — 1) and L x x (L x + 1) 
systems. This gives two results which are presumably 
biased in opposite directions, which can be seen for the 
odd system sizes in Fig. [3] However, our statistical er- 
rors are as large as the bias introduced, so we simply 
average the two results. In comparing sequences of L-2L 
pairs, we check if the result depends on L to extrapolate 
to the thermodynamic limit. We find that the extracted 
values of q c (L, IV) show no visible dependence on sys- 
tem size, as shown in Fig. H) Thus the above estimates 
of q c and v appear to be indicative of the values in the 
thermodynamic limit. In contrast to the two-dimensional 
Ising model, the three-spin model appears to have much 
weaker reentrance, if it is reentrant at all. The phase di- 
agrams for the two models therefore differ considerably 
below the Nishimori line. 



methods are quite limited in scoped and even highly- 
sophisticated heuristic methods can only produce reliable 
results up to L = 14.— 



A. Model 

The two-spin Edwards-Anderson Ising spin glass is 
given by the Hamiltonian 



U 



7 , JijSiSj, 



(10) 



where the spins Si sit at the sites of a cubic lattice and 
pair wise interactions (p = 2) are over nearest neigh- 
bors. Here we consider bimodal disorder, where the bond 
strengths are given by the same probability distribution 
as the plaquettes in Eq. ([2]), 



P{Ja) = (1 - q)8(J l0 - 1) + q5{J l3 + 1). 



(11) 



The zero-temperature fcrromagnet-to-spin-glass tran- 
sition has been investigated for this problem using a so- 
phisticated heuristic algorithm,— where the critical dis- 
order strength was found to be q c = 0.222(5). This work 
did not take into account possible biases in the sampling 
of ground states, which could lead to incorrect results, 
as pointed out in Ref. Eol . The finite-temperature mul- 
ticritical point has been evaluated for this model in a 
number of studies, with the results summarized in Ta- 
ble [III and shown graphically in Fig. [5] Note that, where 
necessary, we have converted estimates so they are all in 
the same terms. In particular, Ref. Hj] quotes only the 
temperature T* of the multicritical point, which is easily 
converted to q* , and the value given for v in Ref. HH is for 
thermal perturbations, with v = l/j/2, while the value we 
quote ) is more relevant to low-temperature studies 
because it describes the disorder perturbations that are 
used at the zero-temperature fixed point. The tempera- 
ture of the multicritical point (MCP) can be determined 
directly from the Nishimori condition 



1 - 2q = tanh(l/T), 



(12) 



with values in the literature in the range T* sa 1.67 — 
1.69. It is believed that the data below the multicritical 
point are universal and in a different universality class 
than at the multicritical point.— This view is consistent 
with the known data. Furthermore, the q c values show 
strong evidence for reentrance in the phase diagram. 



IV. THREE-DIMENSIONAL ISING SPIN GLASS 



B. Results 



As another application of this algorithm, we inves- 
tigate the ferromagnet-to-spin-glass transition in the 
three-dimensional Edwards Anderson Ising spin glass 
with bimodal disorder. This model is important in the 
study of disordered materials.— It is also notoriously dif- 
ficult to study numerically. It is NP-hard^ so exact 



We use the same code for simulating the three-spin 
Ising spin glass to also find ground states of the p = 2 
Ising spin glass in three dimensions. We have computed 
the Binder ratio in the vicinity of the critical disorder 
strength for systems of sizes 4 3 -10 3 . In the Edwards- 
Anderson model, the restricted Binder ratio is identical 
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TABLE II: Selection of different estimates of the critical 
concentration q c and the critical exponent v for the three- 
dimensional bimodal Edwards- Anderson Ising spin glass. The 
estimates are plotted in Fig.[S]and show some variations. Note 
that the values quoted for v at the multicritical point (MCP) 
correspond to disorder perturbations and are not the same v 
one would measure by only perturbing the temperature T.— 
(The roman numerals are used to tag the different estimates 
in the figure). 



Authors 



T 



Ozeki and Nishimori [Ref. 44, I] MCP 0.233(4) 0.51(6) 

Singh [Ref. 41, II] MCP 0.234(2) 0.85(8) 

Hasenbusch et. al. [Ref. 42, III] MCP 0.23180(4) 0.98(5) 

Ceccarelli et. al. [Ref. 43, IV] 1.0 0.2295(2) 0.91(3) 

Ceccarelli et. al. [Ref. 43, V] 0.5 0.2271(2) 0.96(2) 

Hartmann [Ref. L7, VI] 0.0 0.222(5) 1.1(3) 



This study [VII] 



0.0 



0.2253(7) 1.07(7) 




TABLE III: Summary of simulation parameters for the three- 
dimensional Edwards- Anderson model of size L 3 , with pop- 
ulation size N p and averaging over N a different instances of 
disorder. The smallest [largest] value of q simulated is (/mm 

[<?max] in Steps Of (/stop. 



L 




<?stcp 


ymax 


N p 


N s 


4 


0.20 


0.01 


0.25 


1024 


3000 


6 


0.20 


0.01 


0.25 


3456 


3000 


8 


0.21 


0.01 


0.24 


8192 


1700 


10 


0.21 


0.01 


0.24 


16000 


1200 



1 

0.95 

0.9 
0.85 

0.8 
0.75 

0.7 
0.65 

0.6 
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FIG. 6: (Color online) Scaling collapse of the Binder ratio 
g according to Eq. © for the three-dimensional Edwards 
Anderson model. The dashed line is a curve fit to a third- 
order polynomial. Optimal data collapse is obtained for 
q c = 0.2253(7) and v = 1.07(7). 



FIG. 5: (Color online) Graphical representation of different 
estimates of the critical concentration q c and the critical ex- 
ponent v for the three-dimensional Edwards- Anderson model, 
as shown in Table [TT] (tagged via roman numerals) . The cen- 
ter of an ellipse corresponds to the optimal estimate of q c 
and v, while the size of the ellipse corresponds to the error 
bar attached to it. The red ellipses [I - III] are at the multi- 
critical point, which is believed be in a different universality 
class than the low- and zero-temperature case. The blue el- 
lipses [VI, VII] are zero-temperature estimates (the filled el- 
lipse represents the results from this study) and the purple 
ellipses [IV, V] are for finite temperatures below the multi- 
critical point. The fact that the finite-temperature (red) and 
zero-temperature (blue) estimates do not overlap in the hor- 
izontal direction is indicative of a reentrant behavior in the 
phase diagram of the model. 



to the standard definition of the Binder ratio because 
only even moments of the magnetization are used; here 
9 — 9Ri so the scaling relation in Eq. (JS) is again ex- 
pected to hold. Figure [6] shows the scaling collapse for 
this model. The optimal values of the critical parameters 



extracted from this best fit are 

q c = 0.2253(7) v= 1.07(7), (13) 

in agreement with the results of Ref. [T3- The system 
sizes available do not permit the more thorough finite-size 
scaling analysis shown in Fig. 21 although the corrections 
to scaling appear to be quite weak, suggesting that this 
result is robust. A summary of the simulation parameters 
is shown in Table Unl 

V. CONCLUSIONS 

We have demonstrated the applicability of the genetic 
ground-state algorithm presented in Ref. [l4| to sample 
among the many ground states in highly-degenerate NP- 
hard models. Using this technique, spin configurations 
and therefore magnetizations of typical ground-state con- 
figurations may be accessed, allowing the construction of 
parameters such as the Binder ratio without bias, that 
are effective for determining the location of phase tran- 
sitions. 
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Above the Nishimori line, the phase diagram of the 
three-spin Ising model on a triangular lattice closely re- 
sembles that of the square-lattice Ising model with two- 
spin interactions. Although the models are in different 
universality classes, the measured critical values for the 
pure case and on the Nishimori line are the same. The re- 
sults here show that this similarity does not extend below 
the Nishimori line, that is, the models have quite differ- 
ent transition disorder strengths at T = (sec Fig. [4} . 
We show that q c {T = 0) = 0.1072(1). 

Furthermore, using a modest numerical effort we de- 
termine the ferromagnet-to-spin-glass transition for the 
three-dimensional Edwards-Anderson Ising spin glass, 
that is, q c = 0.2253(7), in agreement with previous 
measurements Comparing against estimates of the 
multicritical point, this three-dimensional model pos- 



sesses reentrance in its phase diagram by an amount 
which is similar to the two-dimensional model (similar 
results were obtained recently in Ref. |H) . 
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